---
title: "Mass survey - by country"
author: "Endre Borbáth, Swen Hutter"
date: "20.07.2025"
output:
  html_document:
    toc: yes
---

```{r, echo = F, message = F, include=FALSE}

# R setup

rm(list = ls()) 

library(foreign)
library(dplyr)
library(tidyr)
library(ggplot2)
library(stringr)
library(haven)
library(forcats)
library(lubridate)
library(texreg)
library(DescTools)
library(broom)
library(ordinal)
library(here)
set.seed(1991)

cntry_ext <- c("DE") # adjust with "IT", "PL" depending on which country you'd like to run it for

```


```{r, include = F, echo = F, message = F}

# Data management

pol_controls <- c("contact_org", "contact_org2", "contact_org_orig") #"lrscale", "polinterest", "trust_state"

socdem_controls <- c("age", "age2", "edu_lev2", "gndr", "income", "region")

covid_specific <- c("home_child", "home_care", "infection", "neg_income_chng")

factors <- c("edu_lev2", "gndr", "region", "help_volunteer_past",
             "help_volunteer2", "help_volunteer",
             "support_unknown", "support_no_volunteer", "support_neighbor",
             "support_famfriends", "income", 
             "contact_org", "contact_org2", "contact_org_orig", 
             "org_membership", "home_child", "home_care", "infection")

dependent_vars <- c("help_volunteer", "help_volunteer2", "help_volunteer_past",
                    "support_unknown", "support_no_volunteer", "support_neighbor",
                    "support_famfriends", "vlntr_time_change")

load(here("study2_survey.Rdata"))

```


```{r, include = F, echo = F, message = F, warning=FALSE}

reg_dat <- reg_dat %>% 
  filter(cntry %in% cntry_ext)

```


# Volunteering/ supporting others - reference: the whole sample

## With control variables

```{r, include = F, echo = F, message = F, warning=FALSE}

# Regression models

M1 <- glm(as.formula(paste0("help_volunteer2 ~ contact_org + org_membership + ",
                            paste(covid_specific, collapse = "+"), "+",
                            paste(socdem_controls, collapse = "+"))),
    data = reg_dat,
    weights = wgt_2,
    family = "binomial")

M2 <- glm(as.formula(paste0("support_unknown ~ contact_org + org_membership + ",
                            paste(covid_specific, collapse = "+"), "+",
                            paste(socdem_controls, collapse = "+"))),
    data = reg_dat,
    weights = wgt_2,
    family = "binomial")

M3 <- glm(as.formula(paste0("support_no_volunteer ~ contact_org + org_membership + ",
                            paste(covid_specific, collapse = "+"), "+",
                            paste(socdem_controls, collapse = "+"))),
    data = reg_dat,
    weights = wgt_2,
    family = "binomial")

M4 <- glm(as.formula(paste0("support_neighbor ~ contact_org + org_membership + ",
                            paste(covid_specific, collapse = "+"), "+",
                            paste(socdem_controls, collapse = "+"))),
    data = reg_dat,
    weights = wgt_2,
    family = "binomial")

M5 <- glm(as.formula(paste0("support_famfriends ~ contact_org + org_membership + ",
                            paste(covid_specific, collapse = "+"), "+",
                            paste(socdem_controls, collapse = "+"))),
    data = reg_dat,
    weights = wgt_2,
    family = "binomial")

```

# Table 3, appendix A

```{r, results='asis', include = T, echo = F, message = F, warning=FALSE}

R2s <- c(PseudoR2(M1, "McFadden"), PseudoR2(M2, "McFadden"),
         PseudoR2(M4, "McFadden"), PseudoR2(M5, "McFadden"))

htmlreg(l=list(M1, M2, M4, M5), doctype = FALSE, center = TRUE, siunitx=TRUE, 
       use.packages=FALSE, booktabs = TRUE, float.pos = "h!",
       caption = "Individual-level differences in solidarity action",
       single.row=TRUE, no.margin = TRUE,
       custom.model.names = c(" ", "strangers",
                              "neighbours", "family/ friends"),
       custom.header = list("Volunteering" = 1, "Helping" = 2,
                            "Helping" = 3, "Helping" = 4),
       custom.gof.rows = list("McFadden Sq(R)" =  R2s),
       custom.coef.map = list(#"(Intercept)" = "Intercept",
                              "contact_org2"="Seldom (ref: never)", 
                              "contact_org3"="Sometimes",
                              "contact_org4"="Often",
                              "infection1" = "Infection (self/ in household)",
                              "home_child1" = "Children at home",
                              "home_care1" = "Care work at home",
                              "neg_income_chng" = "Negative income change",
                              "org_membership1" = "Member in civ. soc. org",
                              "age" = "Age",
                              "age2" = "Sq(Age)",
                              "edu_lev21" = "High education",
                              "gndr2" = "Female",
                              "income2" = "Coping",
                              "income3" = "Difficult",
                              "income4" = "Very difficult",
                              "region2" = "Mid. town",
                              "region3" = "Small town",
                              "region4" = "Countryside"),
       groups = list("Contacted by an org." = 1:3,
                     "Crisis specific" = 4:7,
                     "Socio-demographics" = 8:12,
                     "Income (ref: comfortable)" = 13:15,
                     "Region (ref: big city)" = 16:18),
       include.deviance = FALSE,
       include.loglik = FALSE,
       # file=paste0(fig_path, "general_participation.tex"),
      custom.note	= "%stars")

```

Results are based on the first wave, weighted by the socio-demographic weight. These are logit models. In the case of volunteering, I coded never as 0 and everything else as 1. In the helping questions I coded everyone who did not help as zero.

\newpage

Predicted probability of being contacted by an organization

```{r, include = F, echo = F, message = F, warning=FALSE}

new_dat <- data.frame("age" = mean(reg_dat$age),  
                      "age2" = mean(reg_dat$age2),
                      "home_child" = as.factor(0),
                      "home_care" = as.factor(0),
                      "infection" = as.factor(0),
                      "neg_income_chng" = mean(reg_dat$neg_income_chng),
                      "edu_lev2" = as.factor(0), 
                      "gndr" = as.factor(1), 
                      "income" = as.factor(1), 
                      "region" = as.factor(1),
                      "org_membership" = as.factor(0),
                      "contact_org" = as.factor(seq(1, 4, 1)))

results <- NULL

for (i in paste0("M", seq(1, 5, 1))) {

# Calculate the predicted probabilities

# probs <- as.data.frame(predict(!!(i), newdata = new_dat, type = "response", se.fit = TRUE))
probs <- as.data.frame(do.call("predict", list(as.name(i), newdata = new_dat, type = "response", se.fit = TRUE)))
probs$contact_org <- seq.int(nrow(probs)) 
probs$models <- i
probs$residual.scale <- NULL

results <- rbind(results, probs)

}

```


### Figure 4, main text

```{r, include = T, echo = F, message = F, warning=FALSE, fig.width=5, fig.height=5, fig.align="center"}

results <- results %>%
  filter(models!="M3") %>% 
  mutate(type=case_when(models=="M1" ~ "Volunteering",
                        models=="M2" ~ "Helping strangers",
                        models=="M4" ~ "Helping neighbours",
                        models=="M5" ~ "Helping family or friends")) %>% 
  mutate(type=factor(type, levels=c("Volunteering", "Helping strangers", 
                                   "Helping neighbours", "Helping family or friends"))) %>% 
  mutate(pm=fit,
         pu=fit + se.fit * 1.96, # 95% confidence interval
         pl=fit - se.fit * 1.96) %>% # 95% confidence interval 
  mutate(contact_org=as.factor(contact_org))

ggplot(results, aes(x=contact_org, y=pm, group=1)) +
  geom_point() + 
  geom_line(linetype="dashed") +
  geom_errorbar(aes(ymin=pl, ymax=pu), width=.2,
                 position=position_dodge(0.05)) +
  facet_wrap(~ type, ncol = 2) +
  scale_x_discrete(labels=c("1"="Never", "2"="Seldom", "3"="Sometimes",
                            "4"="Often")) +
  theme_linedraw() +
  ylab("Predicted probability") + ylim(0, 1) +
  xlab("Being contacted by an organization")
  # theme(axis.title.x=element_blank())

# ggsave(filename="margins_whole_sample.png",
#        plot=last_plot(),
#        path=fig_path,
#        scale=1.1,
#        width=5,
#        height=5)

```

# Table 4, Appendix A

```{r, results='asis', include = T, echo = F, message = F, warning=FALSE}

R2s <- c(PseudoR2(M1, "McFadden"), PseudoR2(M2, "McFadden"),
         PseudoR2(M3, "McFadden"))

htmlreg(l=list(M1, M2, M3), doctype = FALSE, center = TRUE, siunitx=TRUE, 
       use.packages=FALSE, booktabs = TRUE, float.pos = "h!",
       caption = "Individual-level differences in solidarity action",
       single.row=TRUE, no.margin = TRUE,
       custom.model.names = c(" ", " ",
                              "(without volunteering)"),
       custom.header = list("Volunteering" = 1, "Helping strangers" = 2,
                            "Helping strangers" = 3),
       custom.gof.rows = list("McFadden Sq(R)" =  R2s),
       custom.coef.map = list(#"(Intercept)" = "Intercept",
                              "contact_org2"="Seldom (ref: never)", 
                              "contact_org3"="Sometimes",
                              "contact_org4"="Often",
                              "infection1" = "Infection (self/ in household)",
                              "home_child1" = "Children at home",
                              "home_care1" = "Care work at home",
                              "neg_income_chng" = "Negative income change",
                              "org_membership1" = "Member in civ. soc. org",
                              "age" = "Age",
                              "age2" = "Sq(Age)",
                              "edu_lev21" = "High education",
                              "gndr2" = "Female",
                              "income2" = "Coping",
                              "income3" = "Difficult",
                              "income4" = "Very difficult",
                              "region2" = "Mid. town",
                              "region3" = "Small town",
                              "region4" = "Countryside"),
       groups = list("Contacted by an org." = 1:3,
                     "Crisis specific" = 4:7,
                     "Socio-demographics" = 8:12,
                     "Income (ref: comfortable)" = 13:15,
                     "Region (ref: big city)" = 16:18),
       include.deviance = FALSE,
       include.loglik = FALSE,
       # file=paste0(fig_path, "general_participation.tex"),
      custom.note	= "%stars")

```

Results are based on the first wave, weighted by the socio-demographic weight. These are logit models. In the case of volunteering, I coded never as 0 and everything else as 1. In the helping questions I coded everyone who did not help as zero.

\newpage

Predicted probability of being contacted by an organization

```{r, include = F, echo = F, message = F, warning=FALSE}

new_dat <- data.frame("age" = mean(reg_dat$age),  
                      "age2" = mean(reg_dat$age2),
                      "home_child" = as.factor(0),
                      "home_care" = as.factor(0),
                      "infection" = as.factor(0),
                      "neg_income_chng" = mean(reg_dat$neg_income_chng),
                      "edu_lev2" = as.factor(0), 
                      "gndr" = as.factor(1), 
                      "income" = as.factor(1), 
                      "region" = as.factor(1),
                      "org_membership" = as.factor(0),
                      "contact_org" = as.factor(seq(1, 4, 1)))

results <- NULL

for (i in paste0("M", seq(1, 3, 1))) {

# Calculate the predicted probabilities

# probs <- as.data.frame(predict(!!(i), newdata = new_dat, type = "response", se.fit = TRUE))
probs <- as.data.frame(do.call("predict", list(as.name(i), newdata = new_dat, type = "response", se.fit = TRUE)))
probs$contact_org <- seq.int(nrow(probs)) 
probs$models <- i
probs$residual.scale <- NULL

results <- rbind(results, probs)

}

```

### Figure 3, appendix A

```{r, include = T, echo = F, message = F, warning=FALSE, fig.width=5, fig.height=5, fig.align="center"}

results <- results %>%
  mutate(type=case_when(models=="M1" ~ "Volunteering",
                        models=="M2" ~ "Helping strangers",
                        models=="M3" ~ "Helping strangers (without volunteering)")) %>% 
  mutate(type=factor(type, levels=c("Volunteering", "Helping strangers", 
                                    "Helping strangers (without volunteering)"))) %>% 
  mutate(pm=fit,
         pu=fit + se.fit * 1.96, # 95% confidence interval
         pl=fit - se.fit * 1.96) %>% # 95% confidence interval 
  mutate(contact_org=as.factor(contact_org))

ggplot(results, aes(x=contact_org, y=pm, group=1)) +
  geom_point() + 
  geom_line(linetype="dashed") +
  geom_errorbar(aes(ymin=pl, ymax=pu), width=.2,
                 position=position_dodge(0.05)) +
  facet_wrap(~ type, ncol = 3) +
  scale_x_discrete(labels=c("1"="Never", "2"="Seldom", "3"="Sometimes",
                            "4"="Often")) +
  theme_linedraw() +
  ylab("Predicted probability") + ylim(0, 1) +
  xlab("Being contacted by an organization")
  # theme(axis.title.x=element_blank())

# ggsave(filename="margins_whole_sample_appendix.png",
#        plot=last_plot(),
#        path=fig_path,
#        width=6.5,
#        height=3,
#        scale=1.2,
#        dpi=400)

```


# Being contacted and change in the level of volunteering as DVs

## Being contacted by an organization as a DV

```{r, include = F, echo = F, message = F, warning=FALSE}

M1 <- lm(as.formula(paste0("as.numeric(contact_org_orig) ~ ",
                            paste(covid_specific, collapse = "+"), "+",
                            paste(socdem_controls, collapse = "+"))),
    data = subset(reg_dat, org_membership==0),
    weights = wgt_2)

M2 <- lm(as.formula(paste0("as.numeric(contact_org_orig) ~ ",
                            paste(covid_specific, collapse = "+"), "+",
                            paste(socdem_controls, collapse = "+"))),
    data = subset(reg_dat, org_membership==1),
    weights = wgt_2)

```

## Table 1, Appendix A

```{r, results='asis', include = T, echo = F, message = F, warning=FALSE}

htmlreg(l=list(M1, M2), doctype = FALSE, center = TRUE, siunitx=TRUE, 
       use.packages=FALSE, booktabs = TRUE, float.pos = "h!",
       caption = "Individual-level differences in being contacted",
       single.row=TRUE, no.margin = TRUE,
       custom.model.names = c("Non members", "Members"),
       custom.coef.map = list("infection1" = "Infection (self/ in household)",
                              "home_child1" = "Children at home",
                              "home_care1" = "Care work at home",
                              "neg_income_chng" = "Negative income change",
                              "age" = "Age",
                              "age2" = "Sq(Age)",
                              "edu_lev21" = "High education",
                              "gndr2" = "Female",
                              "income2" = "Coping",
                              "income3" = "Difficult",
                              "income4" = "Very difficult",
                              "region2" = "Mid. town",
                              "region3" = "Small town",
                              "region4" = "Countryside"),
       groups = list("Crisis specific" = 1:4,
                     "Socio-demographics" = 5:8,
                     "Income (ref: comfortable)" = 9:11,
                     "Region (ref: big city)" = 12:14),
       include.deviance = FALSE,
       include.loglik = FALSE,
       # file=paste0(fig_path, "general_participation.tex"),
      custom.note	= "%stars")

```

Results are based on the first wave, weighted by the socio-demographic weight. 
An effect of 1 means a 1 point step up on the five points scale above, it's a different interpretation than in the case of probabilities in the logit model above. 
I include below the effect plot.

\newpage

## Figure 3, main text

```{r, include = F, echo = F, message = F, warning=FALSE}

mod1_res <- tidy(M1)
mod1_res <- mod1_res %>% 
  mutate(type="No member")
mod1_ci <- as.data.frame(confint(M1))
mod1_ci$term <- rownames(mod1_ci)
mod1_res <- merge(mod1_ci, mod1_res, by="term", all=TRUE)
mod1_ci <- as.data.frame(confint(M1, level=0.9))
mod1_ci$term <- rownames(mod1_ci)
mod1_res <- merge(mod1_ci, mod1_res, by="term", all=TRUE)

mod2_res <- tidy(M2)
mod2_res <- mod2_res %>% 
  mutate(type="Member")
mod2_ci <- as.data.frame(confint(M2))
mod2_ci$term <- rownames(mod2_ci)
mod2_res <- merge(mod2_ci, mod2_res, by="term", all=TRUE)
mod2_ci <- as.data.frame(confint(M2, level=0.9))
mod2_ci$term <- rownames(mod2_ci)
mod2_res <- merge(mod2_ci, mod2_res, by="term", all=TRUE)

probs <- bind_rows(mod1_res, mod2_res)

```


```{r, include = T, echo = F, message = F, warning=FALSE, fig.width=5, fig.height=5}

plot_dat <- probs %>% 
  filter(term!="(Intercept)") %>% 
  filter(!grepl("cntry", term)) %>% 
  mutate(term=case_when(term=="edu_lev21" ~ "High educ.",
                        term=="gndr2" ~ "Female",
                        term=="region2" ~ "Mid. town",
                        term=="region3" ~ "Small town",
                        term=="region4" ~ "Countryside",
                        term=="age" ~ "Age",
                        term=="age2" ~ "Sq(Age)",
                        term=="infection1" ~ "Infection (self/ in household)",
                        term=="home_child1" ~ "Children at home",
                        term=="home_care1" ~ "Care work at home",
                        term=="neg_income_chng" ~ "Negative income change",
                        term=="income2" ~ "Coping",
                        term=="income3" ~ "Difficult",
                        term=="income4" ~ "Very difficult",
                        TRUE ~ term)) 

titles <- plot_dat %>%
  filter(term %in% c("Age", "Coping", "Difficult", "Female")) %>%  # I just need some unique values
  mutate_at(vars(c(2:9)), ~ NA)

titles$term <- rep(c("Crisis specific", "Socio-demographics", "Income (ref: comfortable)", "Region (ref: big city)"), 1)

plot_dat <- bind_rows(plot_dat, titles) %>% 
  mutate(term=factor(term, levels = c("Crisis specific",
                                      "Infection (self/ in household)",
                                      "Children at home",
                                      "Care work at home",
                                      "Negative income change",
                                      "Socio-demographics",
                                      "Age", "Sq(Age)", "High educ.", 
                                      "Female", 
                                      "Income (ref: comfortable)",
                                      "Coping", "Difficult", "Very difficult",
                                      "Region (ref: big city)",
                                      "Mid. town", "Small town", "Countryside"))) %>%
  mutate(term=fct_rev(term)) %>% 
  mutate(new_labels = as.character(term)) 

space_between_bars <- 0.5

res1 <- ggplot(plot_dat, aes(y=term, color=type)) +
  geom_vline(xintercept = 0, linetype="dashed", color="gray60") +
  geom_point(aes(x=estimate, shape=type), position=position_dodge(width=space_between_bars), size=2) +
  geom_linerange(aes( xmin=`2.5 %`, xmax=`97.5 %`, 
                     group=type), size=0.5, position=position_dodge(width=space_between_bars), key_glyph = "path") +
  geom_linerange(aes(xmin=`5 %`, xmax=`95 %`,
                     group=type), size=1.1, position=position_dodge(width=space_between_bars), key_glyph = "path") +
  xlab("OLS Estimates") +
  scale_color_manual(name="", values=c("#758fc7", "#cc592d"), guide=guide_legend(reverse = TRUE, nrow=1, byrow=TRUE)) +
  scale_shape_discrete(name="", guide=guide_legend(reverse = TRUE, nrow=1, byrow=TRUE)) +
  scale_y_discrete(labels = function(y) stringr::str_wrap(y, width = 22)) +
  theme_linedraw() +
  theme(legend.title=element_blank(), 
        legend.position="bottom",
        legend.key.width = unit(1,"cm"),
        axis.title.y = element_blank(),
        axis.text.y=element_text(face=c(rep('plain', 3), 'bold', 
                                        rep('plain', 3), 'bold', 
                                        rep('plain', 4), 'bold', 
                                        rep('plain', 4), 'bold'),
                                 hjust=1),
        # axis.text.y = element_text(face = c(rep(2, 4))),
        axis.ticks.y = element_blank())

res1

# ggsave(plot=last_plot(),
#        filename="contact_by_membership.png",
#        path = fig_path,
#        width=5,
#        height=5,
#        scale=3,
#        dpi=400,
#        units="cm")



```

\newpage

I also do now the regression to see how do volunteers and those who supported previously unknown people differ from the general population.. In this part I only work with volunteering and supporting unknown people as the dependent variable - again dichotomized. To get to a sample of the ones that were reached by organizations, I dichotomize the contact variable and distinguish those who report some form of contact (one third of the total sample). I then split the sample and run the same model both on the full sample and on the subset of those who were contacted.

### Table 2, Appendix A

\newpage

## Volunteering

```{r, include = F, echo = F, message = F, warning=FALSE}

# Regression models

M1 <- glm(as.formula(paste0("help_volunteer2 ~ org_membership + ",
                            paste(covid_specific, collapse = "+"), "+",
                            paste(socdem_controls, collapse = "+"))),
    data = subset(reg_dat, as.numeric(contact_org)==1),
    weights = wgt_2,
    family = "binomial")

M2 <- glm(as.formula(paste0("help_volunteer2 ~ org_membership + ",
                            paste(covid_specific, collapse = "+"), "+",
                            paste(socdem_controls, collapse = "+"))),
    data = subset(reg_dat, as.numeric(contact_org)>1),
    weights = wgt_2,
    family = "binomial")

```


```{r, results='asis', include = T, echo = F, message = F, warning=FALSE}

R2s <- c(PseudoR2(M1, "McFadden"), PseudoR2(M2, "McFadden"))

htmlreg(l=list(M1, M2), doctype = FALSE, center = TRUE, siunitx=TRUE, 
       use.packages=FALSE, booktabs = TRUE, float.pos = "h!",
       caption = "Individual-level differences in volunteering",
       single.row=TRUE, no.margin = TRUE,
       custom.model.names = c("Not contacted", "Contacted"),
       custom.gof.rows = list("McFadden Sq(R)" =  R2s),
       custom.coef.map = list("infection1" = "Infection (self/ in household)",
                              "home_child1" = "Children at home",
                              "home_care1" = "Care work at home",
                              "neg_income_chng" = "Negative income change",
                              "org_membership1" = "Member in civ. soc. org",
                              "age" = "Age",
                              "age2" = "Sq(Age)",
                              "edu_lev21" = "High education",
                              "gndr2" = "Female",
                              "income2" = "Coping",
                              "income3" = "Difficult",
                              "income4" = "Very difficult",
                              "region2" = "Mid. town",
                              "region3" = "Small town",
                              "region4" = "Countryside"),
       groups = list("Crisis specific" = 1:4,
                     "Socio-demographics" = 5:9,
                     "Income (ref: comfortable)" = 10:12,
                     "Region (ref: big city)" = 13:15),
       include.deviance = FALSE,
       include.loglik = FALSE,
       # file=paste0(fig_path, "general_participation.tex"),
      custom.note	= "%stars")

```

Results are based on the first wave, weighted by the socio-demographic weight. Let's see the figure with the odds ratios:

\newpage

### Figure 2, Appendix A

```{r, include = F, echo = F, message = F, warning=FALSE}

mod1_res <- tidy(M1)
mod1_res <- mod1_res %>% 
  mutate(type="Not contacted")
mod1_ci <- as.data.frame(confint(M1))
mod1_ci$term <- rownames(mod1_ci)
mod1_res <- merge(mod1_ci, mod1_res, by="term", all=TRUE)
mod1_ci <- as.data.frame(confint(M1, level=0.9))
mod1_ci$term <- rownames(mod1_ci)
mod1_res <- merge(mod1_ci, mod1_res, by="term", all=TRUE)

mod2_res <- tidy(M2)
mod2_res <- mod2_res %>% 
  mutate(type="Contacted")
mod2_ci <- as.data.frame(confint(M2))
mod2_ci$term <- rownames(mod2_ci)
mod2_res <- merge(mod2_ci, mod2_res, by="term", all=TRUE)
mod2_ci <- as.data.frame(confint(M2, level=0.9))
mod2_ci$term <- rownames(mod2_ci)
mod2_res <- merge(mod2_ci, mod2_res, by="term", all=TRUE)

probs <- bind_rows(mod1_res, mod2_res)

```


```{r, include = T, echo = F, message = F, warning=FALSE, fig.width=5, fig.height=5}

plot_dat <- probs %>% 
  filter(term!="(Intercept)") %>% 
  filter(!grepl("cntry", term)) %>% 
  mutate(odds = exp(estimate),
         odds.outter.upper = exp(`2.5 %`),
         odds.outter.lower = exp(`97.5 %`),
         odds.inner.upper = exp(`5 %`),
         odds.inner.lower = exp(`95 %`)) %>% 
  mutate(term=case_when(term=="edu_lev21" ~ "High educ.",
                        term=="gndr2" ~ "Female",
                        term=="region2" ~ "Mid. town",
                        term=="region3" ~ "Small town",
                        term=="region4" ~ "Countryside",
                        term=="org_membership1" ~ "Memb. in civ. soc. org",
                        term=="age" ~ "Age",
                        term=="age2" ~ "Sq(Age)",
                        term=="infection1" ~ "Infection (self/ in household)",
                        term=="home_child1" ~ "Children at home",
                        term=="home_care1" ~ "Care work at home",
                        term=="neg_income_chng" ~ "Negative income change",
                        term=="income2" ~ "Coping",
                        term=="income3" ~ "Difficult",
                        term=="income4" ~ "Very difficult",
                        TRUE ~ term)) 

titles <- plot_dat %>%
  filter(term %in% c("Age", "Coping", "Difficult", "Female")) %>%  # I just need some unique values
  mutate_at(vars(c(2:9), c(11:15)), ~ NA)

titles$term <- rep(c("Crisis specific", "Socio-demographics", "Income (ref: comfortable)", "Region (ref: big city)"), 1)

plot_dat <- bind_rows(plot_dat, titles) %>% 
  mutate(term=factor(term, levels = c("Crisis specific",
                                      "Infection (self/ in household)",
                                      "Children at home",
                                      "Care work at home",
                                      "Negative income change",
                                      "Socio-demographics",
                                      "Memb. in civ. soc. org",
                                      "Age", "Sq(Age)", "High educ.", 
                                      "Female", 
                                      "Income (ref: comfortable)",
                                      "Coping", "Difficult", "Very difficult",
                                      "Region (ref: big city)",
                                      "Mid. town", "Small town", "Countryside"))) %>%
  mutate(term=fct_rev(term)) %>% 
  mutate(new_labels = as.character(term)) %>% 
  mutate(DV_name = "Volunteering") 

space_between_bars <- 0.5

res1 <- ggplot(plot_dat, aes(y=term, color=type)) +
  geom_vline(xintercept = 1, linetype="dashed", color="gray60") +
  geom_point(aes(x=odds, shape=type), position=position_dodge(width=space_between_bars), size=2) +
  geom_linerange(aes(xmin=odds.outter.lower, xmax=odds.outter.upper, 
                     group=type), size=0.5, position=position_dodge(width=space_between_bars), key_glyph = "path") +
  geom_linerange(aes(xmin=odds.inner.lower, xmax=odds.inner.upper, 
                     group=type), size=1.1, position=position_dodge(width=space_between_bars), key_glyph = "path") +
  xlab("Odds Ratio") +
  scale_color_manual(name="", values=c("#758fc7", "#cc592d"), guide=guide_legend(reverse = TRUE, nrow=1, byrow=TRUE)) +
  scale_shape_discrete(name="", guide=guide_legend(reverse = TRUE, nrow=1, byrow=TRUE)) +
  scale_y_discrete(labels = function(y) stringr::str_wrap(y, width = 22)) +
  theme_linedraw() +
  theme(legend.title=element_blank(), 
        legend.position="bottom",
        legend.key.width = unit(1,"cm"),
        axis.title.y = element_blank(),
        axis.text.y=element_text(face=c(rep('plain', 3), 'bold', 
                                        rep('plain', 3), 'bold', 
                                        rep('plain', 5), 'bold', 
                                        rep('plain', 4), 'bold'),
                                 hjust=1),
        # axis.text.y = element_text(face = c(rep(2, 4))),
        axis.ticks.y = element_blank())

# ggsave(plot=last_plot(),
#        filename="volunteering_by_contact.png",
#        path = fig_path,
#        width=5,
#        height=5,
#        scale=3,
#        dpi=400,
#        units="cm")

res1 <- res1 +
  facet_wrap(~DV_name)

res1

```

